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ABSTRACT 


The  vertical  transport  of  heat  by  the  diffusive  layer 
in  the  Arctic  thermocline  is  a  critical  element  of  the  high- 
latitude  climate,  and  yet,  after  decades  of  research,  the 
extant  estimates  remain  highly  controversial.  Laboratory- 
based  estimates  of  vertical  heat  fluxes  originating  from  the 
thermohaline  staircases  of  the  thermocline  are  typically  on 
the  order  of  0.1W/m2  .  This  study  suggests  that  these 
laboratory  experiments  underestimate  the  vertical  heat 
fluxes  and  exceed  their  calculations  by  nearly  an  order  of 
magnitude . 

We  first  guantify  the  typical  density  ratio,  step 
height  and  temperature  gradient  within  the  diffusive 
staircases  of  the  Beaufort  Gyre.  Then,  these  characteristics 
are  used  as  an  input  into  a  numerical  model,  which  simulates 
the  vertical  heat  fluxes  driven  by  the  double  diffusive 
processes.  The  series  of  two-dimensional  simulation  runs 
consistently  calculated  heat  fluxes  on  the  order  of  1 W / m2  . 
In  addition,  analysis  of  a  three-dimensional  simulation 
suggests  that  the  three-dimensional  fluxes  substantially 
exceed  their  two-dimensional  counterparts.  A  detailed 
analysis  of  the  laboratory  measurements  suggests  that  the 
empirical  coefficients  estimated  scaling  factors  from  these 
experiments  are  inconsistent  with  the  corresponding 
numerical  simulations.  These  findings  suggest  that 
laboratory  derived  flux  laws  cannot  be  directly  applied  to 
the  Arctic  Ocean  and  that  further  investigations  into 
double-diffusive  convective  processes  are  warranted. 
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I. 


INTRODUCTION 


In  the  February  2009  issue  of  U.S.  Naval  Institute's 
Proceedings  Magazine,  Rear  Admiral  David  Gove,  U.S.  Navy, 
began  his  article  on  Arctic  Melt  by  stating  "Changes  in  the 
Arctic  environment  —  no  matter  the  cause  —  are  a  great 
national  security  concern."  He  further  asserted: 

To  ensure  complete  maritime  domain  awareness  in 
the  region,  and  to  provide  our  forces  a 
competitive  advantage,  it  will  be  necessary  to 
have  a  comprehensive  knowledge  of  the  physical 
environment ...  Computer-based  ocean  and 
atmospheric  models  must  be  adjusted  to  the 
geophysical  peculiarities  of  high  latitudes. 

In  December  2007,  the  United  States  Environmental 
Protection  Agency  reported  that  warming  rates  in  the  Arctic 
will  be  the  greatest  of  any  other  region  in  the  world,  due 
in  part  to  the  retreating  of  ice  formations.  Less  ice  on 
the  surface  of  the  ocean  increases  the  amount  of  the  sun' s 
energy  that  is  absorbed  into  the  water  column,  thereby 
further  warming  the  planet.  Add  to  this  scenario  the 
possibility  that  the  warmer  Atlantic  Deep  Water  sliding 
beneath  the  colder  and  fresher  Upper  Polar  Deep  Water  (see 
Figure  1)  can  also  convect  enough  energy  to  add  to  the 
demise  of  the  Arctic  sea  ice. 

A  pronounced  feature  of  the  Arctic  thermocline  is 
related  to  the  abundance  of  the  double-diffusive  structures 
in  the  region.  The  ultimate  goal  of  this  thesis  is 
physically  based  parameterization  and  prediction  of  vertical 
double-diffusive  heat  fluxes.  It  is  equally  important  to 
determine  whether  these  fluxes  are  significant  enough  to 
affect  the  Arctic  Climate.  This  study  is  based  on  a 
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combination  of  data  analysis  and  numerical  modeling.  First, 
we  quantify  the  typical  density  ratio,  step  height  and 
temperature  gradient  within  the  diffusive  staircases  of  the 
Beaufort  Gyre.  Then,  the  data  are  used  as  input  into  a 
numerical  simulation,  which  simulates  the  vertical  heat 
fluxes  within  the  gyre  due  to  double  diffusive  layering 
processes . 


Figure  1.  Cross  section  of  Arctic  water  mass.  Cold  and 
relatively  fresh  water  sits  atop  more  warmer  and 
saline  Atlantic  water.  Illustration  by  Jayne  Doucette, 
Woods  Hole  Oceanographic  Institution .  Obtained  from 

WHOI  ITP  website: 

http : / /www.whoi . edu/ oceanus/viewArticle 


Despite  the  critical  role  played  by  the  Arctic  Ocean  in 
the  global  climate,  the  vertical  heat  and  salt  fluxes 
through  the  main  Arctic  thermocline  remain  a  source  of 
greatest  uncertainty.  Do  these  fluxes  have  a  role  in  the 
Arctic  climate?  In  1965,  Turner  derived  his  well-known  4/3 
flux  law  based  on  dimensional  analysis,  which  attempted  to 
quantify  heat  flux,  FT  : 
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aFT=C(Rp )• 


1 

3  4 

■(aAT)i 


(1) 


v  17  y 


where  C  is  a  function  of  R  is  the  background  density 


ratio,  kt  is  the  molecular  diffusivity  of  heat,  g  is  the 
gravitational  acceleration,  v  is  the  kinematic  viscosity, 
a  is  the  coefficient  of  thermal  expansion,  and  AT  is  the 
temperature  variation.  Others  would  follow  (Padman  & 
Dillon,  1987;  Marmarino  &  Caldwell,  1976;  Kelley,  1990)  and 
attempt  to  calibrate  Turner' s  law  with  empirical  data 
garnered  from  laboratory  experiments. 

The  usual  set  up  typically  involves  introducing  a  two- 
layer  stratification  (cold  and  fresh  water  on  top  of  warm 
and  salty)  in  the  experimental  tank  and  monitoring  T  and  S 
variation  in  the  layers.  Conservation  laws  are  then  used  to 
infer  the  temperature  and  salinity  fluxes.  One  of  the 
assumptions  is  that  there  are  no  heat  losses  to  the 
sidewalls.  Another  and  perhaps  even  more  critical  concern 
with  this  setup  is  the  fact  that  the  interaction  between  the 
lower  layer  and  the  tank  bottom  is  not  an  adequate  model  for 
the  layers  in  nature  (Padman  &  Dillon,  1987)  .  The  rigid 
boundaries  pose  a  hindrance  to  the  convective  mixing  and  as 
we  argue  here,  results  in  a  significant  underestimation  of 
the  vertical  transport.  These  attempts  have  been  criticized 
as  simplistic  and  prone  to  large  errors  and  yielded 
estimates  of  the  vertical  heat  flux  in  the  staircase  region 

on  the  order  of  0.02W/m2to  0.22W/m 1  It  has  been  noted 

that  heat  fluxes  less  that  lW/m2  are  not  significant  enough 
to  affect  the  Arctic  Ocean  surface  heat  budget  (Timmermans 
et  al . ,  2008 ) . 
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Numerical  modeling  is  the  simplest  and  most  obvious 
tool  for  evaluating  the  mixing  characteristics  in  the  Arctic 
Ocean,  and  it  will  be  extensively  used  in  our  study, 
supplemented  by  the  analysis  of  relevant  observations. 
Several  attempts  have  already  been  made  (Molemaker  & 
Dijkstra,  1997;  Prikasky,  2007)  to  address  the  problem 
numerically,  but  computational  difficulties  precluded 
simulations  of  the  parameter  regime  relevant  for  the  Arctic 
staircases.  However,  recent  advances  in  high-performance 
computing  make  it  possible  now  to  perform:  i)  first  three- 
dimensional  simulation  of  the  diffusive  convection  in  the 
numerically  accessible  regime  and  ii)  realistic  two- 
dimensional  simulations  in  the  parameter  range  directly 
applicable  to  the  Arctic  staircases.  The  analysis  in  this 
thesis  is  based  on  both  two-dimensional  and  three- 
dimensional  models.  A  series  of  simulation  runs  are 
performed  in  which  we  systematically  vary  the  background 
parameters,  within  the  limits  relevant  for  the  Beaufort  Gyre 
staircases,  record  and  analyze  the  vertical  heat  fluxes  and 
their  dependencies. 

This  thesis  is  set  up  as  follows.  In  Chapter  II,  we 
introduce  the  double-diffusive  processes  and  discuss  their 
dynamics.  In  Chapter  III,  we  investigate  the  Beaufort  Gyre 
and  calculate  the  parameters  characteristic  of  the  Arctic 
location.  Chapter  IV  discusses  in  depth  the  numerical  model 
and  guantifies  heat  fluxes  calculated  in  both  two-dimensions 
and  three-dimensions.  These  fluxes  are  analyzed  and 
compared  with  laboratory  derived  fluxes.  Finally,  our 
conclusions  are  drawn,  and  recommendations  for  future  work 
in  this  field  of  study  are  offered. 
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II.  DOUBLE -DIFFUSIVE  CONVECTION:  BACKGROUND 


A.  DOUBLE -DIFFUSIVE  CONVECTION 

Turner  (1973)  describes  double-diffusive  convection  as 
a  phenomenon  that  occurs  when  there  are  gradients  of  two  or 
more  properties  with  different  molecular  dif fusivities .  The 
double-diffusive  instability  represents  an  efficient 
mechanism  that  releases  potential  energy  stored  in  one  of 
the  density  components,  even  when  the  density  distribution 
is  statically  stable.  With  regard  to  the  ocean,  net 
differences  in  the  diffusivity  of  heat  and  salt  are  of 
interest.  At  20  degrees  Celsius  the  diffusivity  of  heat 
( kt=  1.49x10~7  m2/s)  is  on  the  order  of  100  times  greater 


than  that  of 

salt 

(  ks=  0 . 12 9x1 0~ 7 

m2/s) 

The  disparity  in 

diffusivities 

is 

even  greater  in 

the 

high  latitude  ocean 

(  kt=  1 . 39xl(T7 

m2/s, 

r  ks  =0 . 00 68x1  CT7 

m2/s) 

(Kraus  &  Businger, 

1994)  .  There  are  two  types  of  double-diffusion,  which  are 
known  as  the  salt  finger  regime  and  the  diffusive  convection 
regime . 

1 .  Salt  Finger  Regime 

Salt  fingers  are  usually  observed  in  warm  climates  such 
as  the  tropics  and  subtropics,  where  the  combination  of 
evaporation  and  solar  heating  at  the  ocean  surface  results 
in  a  layer  of  warm,  salty  water  atop  relatively  cooler  and 
fresher  water.  The  physical  mechanism  of  salt  fingering  is 
illustrated  in  Figure  2.  Imagine  displacing  a  parcel  of 
warm,  salty  water  across  the  thin  interface  into  a  region  of 
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colder  and  fresher  water.  Due  to  the  higher  diffusion  rate 
of  heat,  the  parcel  of  water  quickly  loses  its  heat  to  its 
new  surrounding.  At  the  same  time,  it  diffuses  salt  to  the 
less  saline  environment,  but  100  times  slower.  As  a  result, 
the  parcel  of  water  remains  saltier  and  therefore  denser 
than  the  surrounding  fluid  and  continues  to  sink  further. 
This  instability  is  manifested  by  the  appearance  of  narrow 
and  elongated  vertical  filaments,  a.k.a.  salt  fingers. 


Salt  Fingers 


Figure  2.  Salt  Fingers.  Fingers  of  salty  and  cooler  water 
extend  downward  from  the  interface  due  to  the 
difference  in  diffusivity  rates  of  heat  and  salt. 
(After  Ruddick,  1997) 


2 .  Diffusive  Convection  Regime 

While  salt  fingering  has  been  a  subject  of  considerable 
interest  since  its  discovery,  the  diffusive  convection,  on 
the  other  hand,  historically  received  less  attention.  One 
of  the  reasons  is  related  to  its  prevalence  in  high-latitude 
locations,  where  field  programs  are  more  demanding. 
Nevertheless,  clear  signatures  of  the  diffusive  convection 
have  been  observed  and  recorded  (Kelley,  1990;  Padman  & 
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Dillon,  1987)  .  This  study  is  focused  on  the  diffusive 
convection  regime  where  cold  and  fresh  water  overlies 
relatively  warm  and  salty  water. 

To  understand  the  dynamics  of  diffusive  convection, 
imagine  displacing  a  parcel  of  warm  and  salty  water,  as 
illustrated  in  Figure  3.  This  time,  it  is  displaced  upward 
across  the  thin  interface  into  cooler  and  fresher  water.  As 
expected,  the  parcel  of  water  quickly  diffuses  the  heat  it 
carried  with  it  to  its  new  surroundings  and  slowly  releases 
its  salt.  However,  since  the  parcel  of  water  has  now  cooled 
and  has  not  lost  much  of  its  salinity,  it  sinks  across  the 
interface  and  overshoots  its  original  starting  point  due  to 
increased  density  due  to  heat  loss.  As  it  returns  to  warmer 
waters,  it  once  again  gains  heat  and  eventually  will  rise 
above  the  interface  to  once  again  lose  its  heat.  This  over¬ 
stable  process  is  known  as  oscillatory  diffusive  convection. 
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Figure  3.  Oscillatory  regime.  The  fluid  parcel  is  caught 
oscillating  up  and  down  shuttling  heat  upward. 

This  form  of  instability  is,  however,  of  limited 
interest,  since  diffusive  oscillations  are  eventually 
replaced  by  the  quasi-steady  patterns,  characterized  by  the 
layered  distribution  of  temperature  and  salinity  (Prikasky, 
2007)  .  Therefore,  most  commonly  found  in  the  Arctic  is  the 
mode  of  diffusive  layering.  Well-mixed  layers  separated  by 
sharp  diffusive  interfaces  are  responsible  for  stair-like 
appearance  the  temperature  profiles  in  the  central  Arctic 
thermocline.  Such  profiles  are  referred  to  as  thermohaline 
staircases.  Figure  4,  below,  depicts  one  of  the  steps  in 
the  staircase  where  a  layer  of  cold  and  fresh  water  lies 
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over  more  warm  and  salty  water.  A  warm  layer  diffuses  heat 
through  the  thin  interfacial  layer,  immediately  above  it. 
Particles  located  just  above  the  interface  warm  up  and  rise, 
whereas  particles  immediately  below  the  interface  cool  off 
and  sink.  The  resulting  convection  cells  above  and  below 
the  interface  maintain  a  homogeneous  distribution  of  T-S  in 
each  layer. 


A 

z 


Figure  4.  Typical  thermohaline  staircase.  Heat  convects 
upward  through  the  interface  while  the  salt  transport 
is  limited  by  its  lower  molecular  diffusivity. 


Diffusive  staircases  are  observed  where  the  thermocline 
has  a  negative  temperature  gradient.  Heat  is  transferred 
upward  through  a  sequence  of  adjacent  steps  in  the 
staircase.  The  numerical  model  used  in  this  study  makes  use 
of  periodic  boundary  conditions  illustrated  in  Figure  5. 
These  boundary  conditions  offer  a  fairly  realistic 


9 


representation  of  the  mixing  that  occurs  in  nature, 
particularly  in  comparison  to  the  lab  experiments  described 
above . 


Figure  5.  Schematic  diagram  illustrating  the  infinite  series 
of  identical  steps.  In  following  discussion,  the 
temperature  variation  across  each  step  (AT”)  is  defined 
as  a  difference  in  temperature  between  the  centers  of 
two  adjacent  layers  (After  Radko,  2005). 
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Ill .  PRELIMINARY  CALCULATIONS 


A.  ICE-TETHERED  PROFILER 

The  Ice-Tethered  Profiler  (ITP)  was  designed  for 
deployment  on  the  perennial  sea  ice  in  the  polar  oceans.  It 
brings  together  two  main  types  of  oceanographic  sampling 
technologies;  the  mooring  station  and  the  drifter.  The 
marriage  of  these  two  technologies  allows  the  sampling 
apparatus  to  remain  above  the  ocean  surface  in  constant 
contact  with  shore  facilities  while  also  drifting  about  the 
ocean  along  with  the  ice  floe  it  is  moored  to.  The  ITP  has 
three  main  components:  a  surface  subsystem;  a  weighted, 
plastic- j acketed  wire  tether  suspended  from  the  surface 
unit;  and  a  profiler  that  travels  up  and  down  the  jacketed 
wire  via  a  traction  drive  (see  Table  1  further  ITP 
specifications)  .  The  profiler  reports  its  findings  to  the 
surface  unit  using  an  inductive  modem  after  each  one-way 
traverse.  From  the  surface  unit,  this  data  is  then 
transferred  to  shore-based  facilities  using  an  Iridium 
transmitter . 
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Figure  6.  Ice-Tethered  Profiler  schematic.  Illustration 

obtained  from  WHOI  ITP  website: 
http : / /www.whoi . edu/ oceanus/viewArticle 


Size 


Tether  length 
Mass 


Profiling  range 
Endurance 

Surface  unit  temperature 
specification 
Sensors 


Telemetry 

Data  telemetry 
Power 


Surface  unit — 100  cm  in  length.  66  cm  in  diameter 

Underwater  profiling  unit — 123  cm  in  length.  23  cm  in  diameter  (cylindrical  section  diameter  is  15  cm); 

fits  through  25-cm-diameter  ice  hole 
Up  to  800  m 
Surface  unit — 70  kg 
Profiler — 30  kg 
Tether  (800  m) — 250  kg 
Termination  weight — 114  kg 

Approximately  1  500000  m  on  standard  battery  pack 

2.5  to  3  yr  returning  two  750-m  (one  way)  profiles  per  day  (~2000  total) 

Operates  to  — 35°C  for  prototypes,  resumes  full  operation  after  experiencing  temperatures  down 
to  — 48°C 

Sea-Bird  41-CP  CTD 

Dissolved  oxygen,  photosynthetically  available  radiation  (PAR),  chlorophyll  fluorescence,  optical 
backscatter,  and  carbonic  dissolved  organic  matter  (CDOM)  in  development 
Sea-Bird,  Inc.,  inductive  link  from  profiler  to  surface  unit; 

Iridium  link  to  shore  using  dial-up  data  modem 

Typically  50  Kbytes  per  profile  (totaling  100  Mbytes  over  3  yr) 

Lithium  BCX  DD  battery  packs;  3300  Wh  in  surface  package  (after  derating  for  temperature) 

2500  Wh  in  profiler 


Table  1.  Specifications  of  the  ITP  hardware  and  telemetry 

system.  (After  Krishfield  et  al .  ,  2008) 
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Data  received  at  the  WHOI  shore  facility  are  then 
processed  and  made  publicly  available  on  the  ITP  internet 
website  (http :  / /www.whoi .  edu/itp)  as  soon  as  they  are 
received  (Krishfield  et  al . ,  2008). 

There  are  three  general  types  of  data  sets  that  are 
readily  available  for  download  on  the  WHOI  website:  Level  1 
Raw  Data,  Level  2  Real  Time  Data,  and  Level  3  Archive  Data. 
Level  1  Raw  Data  is  simply  raw  sensor  and  engineering  data 
acquired  by  the  profilers.  The  binary  data  are  unpacked  and 
reformatted  into  Matlab,  one  file  per  profile.  Level  2  Real 
Time  Data  incorporates  filtered  location  data  and 
interpolated  to  the  start  times  of  each  profile.  The 
scientific  and  engineering  data  are  then  averaged  in  2-db 
bins.  No  sensor  response  corrections,  calibrations  or 
editing  are  made  in  Level  2  Real  Time  Data,  other  than  the 
internal  calibrations  of  the  CTD  instruments.  Finally,  Level 
3  Archive  Data  represents  the  best  estimates  of  the  ocean 
properties.  The  data  have  had  sensor  response  corrections, 
regional  conductivity  adjustments  based  on  historical  data, 
and  edits  performed.  These  corrections  and  adjustments  are 
not  made  until  the  ITP  has  completed  its  mission  (Krishfield 
et  al . ,  2008 )  . 

B.  DATA  ANALYSIS 

Data  sets  from  ITPs  1,  2,  3  and  5  were  chosen  based  on 
the  tracks  they  ended  up  traveling  over  while  they  were 
deployed.  Figure  7  is  a  graphic  taken  from  the  WHOI  ITP 
website  and  it  depicts  the  area  the  four  ITPs  covered. 
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Figure  7.  ITP  tracks  through  the  Beaufort  Gyre.  Illustration 

obtained  from  WH0I  ITP  website: 
http : / / www. who i . edu/page . do?pid=23099 


Matlab  routines  were  written  to  process  the  data.  First, 
temperature  (Figure  8)  and  salinity  profiles  (Figure  9)  were 
plotted  for  each  ITP.  To  identify  the  range  of  depths  at 
which  the  preponderance  of  thermohaline  steps  existed, 
temperature  versus  depth  was  plotted  in  groups  of  20-50 
profiles  per  plot  as  shown  in  Figure  10.  Next,  layer 
interfaces  were  identified  utilizing  Matlab  script  which 
looked  for  a  maximum  change  in  temperature.  Once  the 
interfaces  were  identified,  the  script  then  identified  the 
corresponding  temperature,  salinity  and  depth  with  each 

interface  and  marked  the  midpoints  of  each  step  shown  in 
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Figure  11.  These  quantities  were  then  used  to  calculate  the 
background  density  ratio,  R  ,  using  equation  (5)  for  each 
profile  recorded: 


,  /?A S 

p  a  AT 

r 


(5) 


where  /?  is  the  coefficient  of  saline  contraction,  a  is  the 
coefficient  of  thermal  expansion,  and  AS  and  AT  are  the 
variations  in  salinity  and  temperature  across  each  step, 
respectively.  Both  (3  and  a  were  calculated  using  the 
Matlab  seawater  toolbox. 
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Figure  8.  Temperature  profile  for  ITP2 .  244  profiles 

depicted . 
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ITP  2  Salinity  vs  Depth 


Figure  9.  Salinity  profile  for  ITP  2.  244  profiles 

depicted . 


Figure  10.  Typical  thermohaline  staircase  structure  for  ITP2 . 
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Figure  11.  Individual  profile  within  ITP  2  data  set.  Marker 
indicates  the  midpoint  of  each  step  as  identified  by 
routine  written  in  Matlab. 


Table  2  summarizes  the  parameters  that  quantify  the 
thermohaline  staircases  within  the  Beaufort  Gyre.  Of  most 
importance  is  the  calculation  of  Density  Ratio.  This 
quantity  is  a  measure  of  intensity  of  diffusive  convection 
and  is  essential  for  the  correct  estimate  of  heat  fluxes 
using  a  numerical  model. 


Profiler 

Mean  Density 
Ratio 

Mean  Step 
Height, 
meters 

Mean  Temperature 
Gradient 

ITP  1 

4.87  +/-  1.72 

2.70  +/-  1.02 

0.0164  +/-  0.0054 

ITP  2 

4.81  +/-  1.52 

2.54  +/-  0.80 

0.0165  +/-  0.0068 

ITP  3 

5.43  +/-  1.70 

2.42  +/-  0.77 

0.0160  +/-  0.0102 

ITP  5 

4.59  +/-  1.33 

2.46  +/-  0.81 

0.0158  +/-  0.0053 

Table  2.  Mean  density,  step  height,  and  temperature 

gradient  within  the  Beaufort  Gyre. 
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IV.  DIRECT  NUMERICAL  SIMULATIONS 


A.  MODEL  DESCRIPTION 

Direct  Numerical  Simulations  was  conducted  in  both  two 
and  three  dimensions  on  Naval  Postgraduate  School  (NPS) 
high-performance  clusters,  Anastasia  and  Kinglear,  in 
addition  to  Department  of  Defense  (DoD)  supercomputer 
clusters,  Babbage,  Davinci,  and  Midnight.  Low  resolution 
model  runs  were  set  up  on  the  NPS  clusters  dedicating  only 
one  processor  per  run.  High-resolution  runs  committed  16 
processors  per  run  across  the  three  DoD  supercomputer 
clusters  and  NPS's  Kinglear.  Because  of  a  wide  range  of 
spatial  and  temporal  scales,  numerical  modeling  of  diffusive 
staircases  represented  a  formidable  computational  challenge, 
particularly  in  three-dimensions,  and  required  use  of 
parallel  architecture.  The  numerical  model  is  a  fully  de- 
aliased  pseudo-spectral  FORTRAN  code,  described  in  Radko 
(2003),  using  the  Message  Passing  Interface  (MPI)  for  three 
dimensions.  The  numerical  setup  was  based  on  Ice-Tethered 
Profiler  observations  (from  the  Beaufort  Gyre  Exploration 
Program  led  by  Woods  Hole  Oceanographic  Institute)  and  the 
results  were  compared  with  existing  models.  The  two- 

dimensional  study  examined  the  diffusive  convection  in  the 
horizontal  (x)  and  vertical  (z) . 

Following  Radko  and  Stern  (1999),  temperature  and 
salinity  fields  are  decomposed  into  the  linear  basic  state 
( T,S )  and  perturbations  ( T,S )  .  The  Boussinesq  equations  of 

_  J_ 

motion  are  non-dimensionalized  using  d  =  (pv/gat)4  as  the 
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length  scale. 


as  the  time 


v  =  — —  the  velocity  scale,  t=  — 
d  kt 

scale,  where  kt  is  the  molecular  diffusivity  of  heat,  v  is 

the  kinematic  viscosity,  g  is  the  gravitational 

acceleration,  a  is  the  coefficient  of  thermal  expansion,  Tz 

is  the  temperature  gradient  and  aT_d  is  the  scale  for  both 

temperature  and  salinity  perturbations  resulting  in 


—  —  =  -Vp  +  V2v  +  (T-  S)k, 
Pr  dt 


V  •  v  =  0, 


dT 


- w  = 


dt 

dS_ 

dt 


RP  w 


v2t, 

=  z-V2S, 


(2) 


where  Pr  =  v'//rr=13  is  the  Prandtl  number,  z=ks/kt  is  the 
molecular  diffusivity  ratio  (Lewis  number) ,  R  is  the 

background  density  ratio  for  the  basic  uniform  S-T  gradient, 
while  v,  T,  and  S  are  non-dimensional  velocity,  temperature, 
and  salinity. 

The  non-dimensional  and  dimensional  temperature  fluxes 
are  related  by 

FT<iim=FTnon-iim^T  I  dz  (3) 

and  the  corresponding  heat  flux  is 

FH&m  ~  FT dim  F pP  (4) 

where  C  is  the  specific  heat  for  sea  water  and  p  is  the 
density  of  sea  water.  Substituting  (3)  into  (4)  yields  the 
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conversion  factor  between  the  dimensional  heat  flux  and  non- 
dimensional  temperature  flux  of  .0094  W/m1  . 

The  strategy  was  to  focus  on  two-dimensional 
simulations  and  compare  the  results  to  computationally 
accessible  three-dimensional  simulations.  The  ratio  of 
three-dimensional  heat  fluxes  to  two  dimensional  heat  fluxes 
would  then  be  used  to  calibrate  all  two-dimensional 
simulations . 

B.  TWO-DIMENSIONAL  SIMULATION 

1.  General  Characteristics 

In  order  to  gain  an  understanding  of  the  mechanism  of 
diffusive  layering  and  its  dependencies  on  the  governing 
parameters,  two-dimensional  simulations  calculated  heat 
fluxes  as  functions  of  background  density  ratio,  ranging 
from  3  to  6,  and  Lewis  numbers  with  ratios  1/200,  1/100, 
1/50  and  1/25  creating  a  4x4  matrix  data  set  shown  in  Table 
2.  The  simulation  runs  were  initiated  from  a  random 
computer-generated  distribution  of  T  and  S.  The  grid  size 
used  in  the  initial  16  two  dimensional  runs  was  2048x4096. 

A  typical  evolutionary  pattern  observed  in  the  two- 
dimensional  simulations  are  shown  in  Figure  12.  The  two 
layered  simulation  is  initiated  at  rest.  Random  noise  gives 
the  system  a  kick  and  initiates  the  diffusive  process. 
Notice  the  huge  spike  in  the  very  beginning  of  the 
simulation  indicating  the  quick  release  of  potential  energy 
by  the  temperature  field.  As  time  increases  the 
fluctuations  in  the  heat  fluxes  eventually  reach 
equilibrium.  Mean  heat  fluxes  were  determined  from  the 
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equilibrated  portion  of  the  graph  (typically  achieved  during 
last  one-third  of  time  record) . 

Figures  13,  14,  15  and  16  are  successive  full 
temperature  field  plots  revealing  the  spatial  structure  of 
the  diffusive  convection.  The  two  layered  fluid  is  at  rest 
initially,  as  shown  in  Figure  13.  Figure  14  depicts  signs 
of  initial  convective  activity  as  the  heat  exchange  begins 
to  diffuse  between  the  two  layers  at  t=50  minutes. 


Figure  12.  Typical  two-dimensional  plot  of  heat  flux  versus 

time . 
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Figure  13. 

grid 


Two-dimensional  temperature  field, 
size  of  2048x4096  at  t=0  minutes. 
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Figure  14.  Two-dimensional  temperature  field.  Plot  is  on  XZ 
grid  size  of  2048x4096  at  t=50  minutes.  Convection 
plumes  become  visible  in  both  layers.  High  values  of 
T  are  shown  in  red.  Low  values  in  blue. 
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The  heat  signatures  of  active  heat  transport  are 
clearly  visible  at  t=200  minutes  in  Figure  15.  This  figure 
represents  an  enlarged  view  of  the  two  fluid  interface  to 
show  more  detail  of  the  dynamic  event. 


Figure  15.  Two-dimensional  temperature  field  at  t=200 

minutes.  View  is  zoomed  in  on  the  interface  to  show 
more  detail  of  convective  dynamic.  Temperature  field 
is  well  mixed  above  and  below  the  interface. 
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Figure  16.  Two-dimensional  temperature  field  plot  at  t=400 
minutes.  View  is  zoomed  in  on  the  interface  to  show 
more  detail  of  convective  dynamic.  Interface  is 
beginning  to  show  signs  of  further  convective  activity 
with  small  plumes  emanating  from  the  right  side  of  the 

picture . 


2 .  Dependence  of  Fluxes  on  the  Density  Ratio 

As  indicated  by  equation  (1),  Turner  suggests  that  heat 
fluxes  are  dependent  on  the  density  ratio  for  a  fluid  with  a 
given  temperature  variation  across  the  interface  as  deduced 
from  the  foregoing  numerical  experiments.  Figure  17 
illustrates  this  dependence  on  density  ratio.  As  expected, 
heat  fluxes  decreased  as  the  density  ratio  increased.  In 
fact,  heat  fluxes  were  doubled  when  density  ratio  increased 
from  Rp=  3  to  Rp=  6. 
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Figure  17.  Two-dimensional  heat  flux  versus  density  ratio. 

Heat  fluxes  decrease  as  density  ratio  increases. 

3 .  Dependence  of  Fluxes  on  the  Lewis  Number 

The  Lewis  number  is  a  non-dimensional  ratio  between  the 
dif fusivities  between  salt  and  heat.  A  study  conducted  by 
Takao  and  Narusawa  (1980)  quantified  the  relationship  of 
Lewis  number  and  heat  fluxes  based  on  a  series  of  laboratory 
experiments  with  substances  of  different  molecular 
diffusivity  ratios.  They  found  that  heat  fluxes  increased 
as  the  Lewis  number  decreased.  This  dependence  on  the  Lewis 
number,  although  less  pronounced,  is  also  found  in  this 
study.  Heat  fluxes  increased  as  Lewis  number  decreased  over 
a  range  from  r=l/25  to  r=l/200  as  depicted  in  Figure  18. 
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Figure  18.  Two-dimensional  heat  flux  versus  Lewis  number. 

Heat  fluxes  decrease  as  Lewis  number  increases. 

The  two-dimensional  simulation  results  are  summarized 
in  Table  3.  The  results  show  that  diffusive  heat  fluxes 
increase  by  a  factor  of  two  to  three  as  the  Lewis  number 
decreases  from  1/25  to  1/200.  It  should  be  noted  that  Lewis 
numbers  on  the  order  of  1/100  to  1/200  are  typical  for  the 
ocean.  As  stated  earlier  in  the  paper,  there  is  a  general 
consensus  that  heat  fluxes  that  are  considerably  less  than 
1  W/m2  are  not  significant  enough  to  affect  the  Arctic  Ocean 
surface  heat  budget.  Here,  the  simulation  consistently 
achieved  heat  fluxes  on  the  order  of  lW/m2  for 
aforementioned  density  ratio  and  Lewis  number  ranges. 
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T 

3 

4 

5 

6 

1/200 

i .  i 

0.78 

0.62 

0.51 

W/m 2 

W/m2 

W/m2 

W/m2 

1/100 

0.91 

0.80 

0.49 

0.40 

W/m 2 

W/m2 

W/m 2 

W/m2 

1/50 

0.71 

0.58 

0.39 

0.25 

W/m1 

W/m 2 

W/m 2 

W/m 2 

1/25 

0.41 

0.30 

0.23 

0.19 

W/m2 

W/m2 

W/m2 

W/m2 

Table  3.  Two-dimensional  diffusive  heat  flux  as  a  function 

of  density  ratio  and  Lewis  number. 

C.  THREE-DIMENSIONAL  SIMULATION 

In  order  to  comply  with  the  limits  set  by  the  available 
computing  resources,  a  three-dimensional  simulation  was  run 
with  R  =  6  and  r=l/10  was  conducted  and  compared  with  the 

corresponding  two-dimensional  run.  This  Lewis  number, 
r=l/10,  was  chosen  due  to  it  being  computationally 
accessible  in  the  three-dimensional  simulation.  Choosing 
smaller  Lewis  numbers  necessitated  resolving  the  salt 
dif fusivities  at  very  fine  scales,  which  proved  to  be 
computationally  prohibitive.  A  grid  size  of  256x512  and 
256x256x512  were  used  for  the  additional  two-dimensional  and 
comparable  three-dimensional  runs,  respectively. 

Figure  19  illustrates  the  simulation  conducted  in 
three-dimensions.  Again,  the  simulation  begins  at  rest  and 
is  perturbed  by  small  amplitude  random  noise.  The  system 
had  an  initial  spike  and  began  to  settle  down  as  time 
increased . 
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Figure  19.  Three-dimensional  heat  flux  versus  time.  R  =  6 

and  r=l/10 


Figures  20  to  23  illustrate  the  evolutionary  process  of 
the  temperature  field  in  three-dimensions.  Notice  the 
"worm-like"  patterns  in  the  horizontal  plane  of  Figure  21 
indicating  regions  of  convective  plumes.  These  patterns 
have  grown  into  larger  strands  of  relatively  cooler 
temperatures  surrounded  by  pools  of  higher  temperatures 
indicating  areas  of  high  and  low  heat  flux,  respectively. 
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1=1/10,  R'=6 


Figure  20.  Three-dimensional  temperature  field  plotted  on  XYZ 
grid  size  of  256x256x512  at  t=0  minutes.  High  values 
of  T  are  shown  in  red.  Low  values  in  blue. 


'  t=1/10,  R'=6 

9 


Figure  21.  Three-dimensional  temperature  field  plotted  on  XYZ 
grid  size  of  256x256x512  at  t=70  minutes.  Convection 
plumes  become  visible  on  vertical  panels  and  are 
indicated  by  lighter  "worm-like"  patterns  in  the 
horizontal  plane.  High  values  of  T  are  shown  in  red. 

Low  values  in  blue. 
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t=1/10,  R  =6 


Figure  22.  Three-dimensional  temperature  field  plotted  on  XYZ 
grid  size  of  256x256x512  at  t=500  minutes.  Plumes  are 
beginning  to  dissipate  into  mixing  region  above  and 
below  the  layer.  High  values  of  T  are  shown  in  red. 

Low  values  in  blue. 


t=1/10,  R  =6 

9 


Figure  23.  Three-dimensional  temperature  field  plotted  on  XYZ 
grid  size  of  256x256x512  at  t=700  minutes.  The 
temperature  field  is  now  well  mixed  above  and  below 
the  interface.  High  values  of  T  are  shown  in  red.  Low 

values  in  blue. 
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D .  CALIBRATION  PROCEDURES 

Resolving  scale  of  finer  detail  becomes  the  issue  when 
attempting  to  realistically  model  the  ocean  in  three- 
dimensions.  Lewis  numbers  on  the  order  of  r=l/200  to 
r=l/50  proved  to  be  computationally  difficult.  However, 
this  range  of  Lewis  numbers  is  more  accessible  when 
simulated  in  two-dimensions.  Figure  24  illustrates  a  2- 
dimensional  plot  of  heat  fluxes  utilizing  the  same  input 
values,  R  =  6  and  r=l/10,  used  in  the  three-dimensional 

run . 


Figure  24.  Two-dimensional  heat  flux  versus  time.  R  =  6  and 

r =1/10 . 
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Figure  25  presents  a  superposition  of  two-dimensional 
and  three-dimensional  simulations  in  Figures  23  and  24  and 
shows  how  heat  fluxes  simulated  in  two-dimensions  compares 
to  those  calculated  in  three-dimensions.  The  mean  value  for 

the  two-dimensional  case  was  0.1  Wjm2  while  the  three- 

r 

dimensional  case  was  larger  by  50%.  Given  this  ratio  of  1.5 
and  assuming  it  is  valid  for  lower  values  of  Lewis  numbers 
in  two-dimensions.  Table  4  indicates  the  calibrated  values 
for  density  ratios  between  3  and  6  with  a  Lewis  number  of 
1/200. 


Figure  25.  Comparison  of  two-dimensional  heat  flux  with 

three-dimensional  heat  flux. 
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*p 

T 

3 

4 

5 

6 

1/200 

i .  i 

0.78 

0.62 

0.51 

W/m 2 

W jm2 

W/m2 

W/m 2 

1/200 

1 . 65 

1 . 17 

0.93 

0.77 

Calibrated 

Wjm2 

Wjm 2 

W/m 2 

W/m 2 

Table  4.  Two-dimensional  diffusive  heat  flux  as  a  function 

of  density  ratio  and  Lewis  number  1/200.  Calibrated 
values  are  50%  larger  than  the  original  based  on  the 
comparison  of  two-dimensional  and  three-dimensional 

simulations . 

E.  COMPARISON  WITH  LABORATORY  EXPERIMENTS 

In  order  to  more  fully  explain  the  difference  in  fluxes 
using  numerical  simulations  and  the  earlier  estimates  in 
Kelley  (1990)  derived  from  laboratory  experiments,  we  now 
examine  the  lab-based  flux  laws  directly.  In  particular,  we 
turn  our  attention  to  the  non-dimensional  coefficient  C  of 
Turner's  4/3  flux  law,  equation  (1)  .  In  Figure  26,  we 

present  the  estimates  of  the  non-dimensional  factor  C  based 
on  the  numerical  experiments  in  Table  4,  superimposed  on  the 
earlier  laboratory  estimates.  Numerical  simulations 

produced  values  of  C  that  were  two  to  ten  times  larger  than 
those  achieved  in  laboratories.  Again,  this  should  raise 
the  question  as  to  the  relevance  of  the  lab  based  flux  laws 
in  application  to  the  Arctic  staircases. 
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Figure  26.  Non-dimensional  coefficient  of  the  4/3  flux  law 
C(R )  from  equation  (1)  .  Data  sources  are  indicated  as 

follows:  dashed  Kelley  (1990),  solid  Marmarino  & 
Caldwell  (1976),  circles  two-dimensional  simulation, 
and  triangles  two-dimensional  values  scaled  by  three- 
dimensional  simulation. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 


The  double-diffusive  estimates  of  the  heat  transport  by 
convection  in  the  Arctic  are  highly  controversial. 
Laboratory-based  estimates  suggest  that  the  vertical  heat 
fluxes  through  the  thermohaline  staircases  of  the 
thermocline  are  typically  on  the  order  of 
0.1  W / m2  (Timmermans,  2008;  Padman  &  Dillon,  1987).  This 
study  indicates  that  these  laboratory  experiments 
underestimate  the  vertical  heat  fluxes  and  the  numerical 
simulations  exceed  laboratory  counterparts  by  an  order  of 
magnitude.  Why  such  a  large  difference  in  values?  It  is  my 
opinion  that  the  laboratory  experiments  do  not  correctly 
model  the  ocean  environment  and  therefore  cannot  accurately 
predict  their  vertical  heat  fluxes.  Laboratory  experiments 
are  conducted  in  small  tanks  that  have  rigid  boundaries, 
quite  different  parameter  ranges  and  are  very  susceptible  to 
experimental  errors.  Numerical  modeling  avoids  many  of 
issues,  allowing  vertical  fluxes  to  pass  through  periodical 
boundaries,  providing  a  more  realistic  model  of  the  oceanic 
conditions.  Heat  fluxes  in  this  study  were  found  to  be  on 
the  order  of  lW/m2  indicating  that  heat  fluxes  that  arise 
from  the  Arctic  thermocline  are  indeed  significant  and  may 
play  a  role  in  the  Arctic  heat  budget. 

The  analysis  in  this  thesis  warrants  future  work  on 
this  subject.  Increasing  the  resolution  of  the  model  can 
imitate  nature  even  more  and  provide  more  realistic  results. 
Engaging  in  the  inverse  modeling  of  the  Arctic  staircases 
will  resolve  the  controversy  caused  by  the  disagreement  of 
different  approaches.  A  further  study  to  develop 
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quantitative  analytical  models  of  diffusive  layering  will 
help  us  better  understand  the  nature  of  the  diffusive 
convection.  Finally,  push  the  limits  of  modeling  in  three- 
dimensions  to  include  lower  Lewis  numbers  and  adding  in  a 
horizontal  T-S  gradient  to,  again,  better  simulate  nature. 

It  is  my  hope  that  these  findings  will  cause  the 
scientific  community  to  re-evaluate  their  methods  and 
simulate  further  research  in  this  field  of  study.  Only  then 
can  the  challenge  set  by  Admiral  Gove  be  met,  improving  our 
prediction  models  and  thereby  improving  our  military' s 
ability  to  improve  our  national  security. 
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